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Abstract 

In this note we study the convergence of monotone P 1 finite element methods on unstructured meshes for fully non- 
linear Harrdlton-Jacobi-Bellman equations arising from stochastic optimal control problems with possibly degenerate, 
isotropic diffusions. Using elliptic projection operators we treat discretisations which violate the consistency conditions 
of the framework by Barles and Souganidis. We obtain strong uniform convergence of the numerical solutions and, under 
non-degeneracy assumptions, strong!^ convergence of the gradients. 



I Introduction 

Hamilton-Jacobi-Bellman (HJB) equations, which are of the form 

-5fy-i-sup(L"i;-d") = 0, (1) 

a 

where the L" are linear first- or second order operators and d" £ L^, characterise the value function of optimal control 
problems. Indeed, one possibility to introduce the notion of solution of (1) is via the underlying optimal control structure. 
An alternative approach is to use the monotonicity properties of the operator which leads to the concept of viscosity 
solutions. While these perceptions are essentially equivalent [17, p.72] both views have been instructive for the design 
and analysis of numerical methods. 

The former approach, based on the discretisation of the optimal control problem before employing the Dynamic Pro- 
gramming Principle, has been proposed in the setting of finite elements in [26, 7, 8], see also the review article [22[ and 
the references therein. Regarding finite difference methods we refer to the book [23] . The latter approach, which is also 
adopted in this note, was firmly established with the contribution [3] by Barles and Souganidis in 1991, providing an ab- 
stract framework for the convergence to viscosity solutions. Starting with [20, 21] techniques were developed to quantify 
the rate of convergence; more recent works are [1, 13]. A third direction was opened by the method of vanishing moments 
which neither enforces discrete maximum principles nor makes use of the underlying optimal control structure but relies 
on a higher order regularisation [16]. For a more comprehensive review of the state-of-the-art in the numerical solution 
of fully non-linear second order equations we refer to [15] . 

In the traditional finite element analysis the multiplicative testing with hat functions is viewed as the discrete analogue 
of the multiplicative testing procedure to define weak solutions of the (variational) differential equation. While elements 
of this viewpoint are implicitly used in Section VII on gradient convergence, we would like to stress a second interpre- 
tation: multiplication with hat functions as regularisation of the residual. Consider for a moment the linear problem 
-a(x)Au{x] = f(x) with smooth functions a and u as well as a hat function <p at the node y^. Let P be the orthogonal pro- 
jection onto the approximation space with respect to the scalar product {v,w) - / Vf ■ Vw dx (given suitable boundary 
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conditions). If y is near then on a fine mesh 

-a(y)AM(j/) = - ja[y)Au[y)(pix)dx~-a(y^) j \uix)<pix)dx- a(y^) j Vu(x)-W^{x)dx- aiy^) j WPuM-W<p{x)dx. 

since := (/>/ ll</>llii(n) approximates a Dirac Delta as the element size is decreased. In contrast, on general meshes, 

-a(y) AM(y) ^ aiy^) j VJm(x) • V<p(.x) dx, nodal interpolant) 

even in the limit as the mesh is refined (see Example 1 below). This indicates that the orthogonality properties of the pro- 
jection of the exact solution into the approximation space play an important role for the understanding of the (pointwise) 
consistency of the finite element scheme. Furthermore, this interpretation may serve as a starting point in selecting a 

discretisation of the HJB operator. 

Our analysis combines the following key elements in a single finite element framework: 

Treatment ofnodally inconsistent discretisations and uniform convergence: The consistency condition (see [3, eqn.(2.4)] 
or [17, p.332]) of Barles and Souganidis is based on a limit involving pointwise values of smooth test functions. This con- 
dition is not satisfied by finite element methods, even for linear equations. Based on an alternative consistency condition 
we show the uniform convergence of finite element solutions to the viscosity solution. 

Gradient convergence: We demonstrate how the coercivity of the linear operators under the supremum is recovered by 
the finite element method in order to control the gradient of the numerical solutions. In a uniformly elliptic setting, this 
leads to strong convergence in ( [0, T] , f/^O) ) . 

Operators of non- negative characteristic form: The presented analysis includes the treatment of partially and fully deter- 
ministic optimal control problems, corresponding to degenerate elliptic operators under the supremum of the Hamilto- 
nian. 

Unstructured meshes: In the spirit of finite element methods the computational domain may be triangulated with an 
unstructured mesh, allowing to capture complex domains more easily than in a finite difference setting. Typically, weaker 
conditions on the mesh than quasi-uniformity can be made. 

Regularisation with second order operators: We highlight that the regularisation with second order elliptic operators is 
sufficient to achieve convergence to the viscosity solution. Indeed, in the example of the method of artificial diffusion, we 
illustrate how the regularisation in the second order fully non-linear case is of the same kind and order as for first order 
linear operators. 

Unconditional time step: Our analysis permits explicit, semi-implicit and fully implicit discretisations in time. FuUy im- 
plicit discretisations in time lead to unconditionally stable schemes. 

The structure of the article is as follows: In Section II we introduce a framework of finite element methods. In Section 
III we study the well-posedness of the discrete systems of equations and describe how these systems are solvable by a 
known globally convergent, locally superlinearly convergent algorithm. Section IV establishes consistency properties 
of elliptic projection operators. This enables us to demonstrate in Section V that the upper and lower envelopes of the 
numerical solutions are sub- and supersolutions. Uniform convergence to the viscosity solution is derived in Section VI 
and is then built upon to analyse the convergence of the gradient in Section VII. We provide a concrete specimen of a 
scheme belonging to our framework by describing the method of artificial diffusion in Section VIII. 

II Problem statement and definition of tlie numerical method 

Let Q be a bounded Lipschitz domain in IR'^, d>2. Let ^ be a compact metric space and 

A — C(n) X C(n,IR'*) X C(n) x C(n), a — {a" ,b" ,c" ,d") 

be continuous, such that the families of functions {a"}aeA. {^'"laeyi. {c"'}a£A and {d"}a^A are equi-continuous. Consider 
the bounded linear operators of non-negative characteristic form [25] 

i": hI{D]^H-^{D], w^-a"Aw + b"-SIw + c"w 
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where a belongs to A. Furthermore, suppose that pointwise d" > 0. Then 

sup||(a",f7",c",d")||c(n)><c(n,Kd)><C(n)><C(n) <oO' ^^P\\L"\\cHa)^c(n) <°°- ^2) 

We assume that the final-time boundary data Vj £ C(n) is non-negative: i^f > on n. For smooth w let 

Hw:=sup[L"w-d"), 

a 

where the supremum is applied pointwise. The HJB equation considered is 

-dtV + Hv = in(0,r)xn, (3a) 
v = on (0, T] X an, (3b) 
v=Vt on{r}xn. (3c) 

Definition 1 ([2, 17]). An upper semi-continuous (respectively lower semi-continuous) function i/ : [0, T] x Q — |R /s a vis- 
cosity subsolution (respectively supersolution) of 

-dtV + Hv = (4) 

on(0, r)xn ifforanyw e. C°°(IRxlR'*) such that v-w has a strict local maximum (respectively minimum) at it, x) e (0, r)xn 

with v{t,x] = w{t,x), gives -dtw{t,x] + Hw[t,x] < 0, (respectively greater than or equal to 0). IfveCiRxW^) of equation 
(4) is a viscosity subsolution and supersolution, then v is called a viscosity solution. 

Let Vi be a sequence of piecewise linear shape-regular finite element spaces with nodes and associated hat fiinctions 
(p^^. Here ( is the index ranging over the nodes of the finite element mesh. Let 1/° be the subspace of functions which 
satisfy homogeneous Dirichlet conditions. It is convenient to assume that e Q for ^ < JV := dim V^; i.e. the index ( first 
ranges over internal and then over external nodes. Set := <^f /||</>f llLi(n)- The mesh size, i.e. the largest diameter of an 
element, is denoted (A;c),-. It is assumed that {^x)i — as i — oo. 

Let hihe the (uniform) time step used in conjunction with V/ , with ^ £ N, and let sj^ be the fcth time step at the refinement 
level i. The set of time steps is 

S,:={5f:fc = 0,...,;|:}. 

Let the ^th entry of d,- w{s^, •) be 

idi u;(5f , ■))( = ' 

hi 

For each a and i find an approximate splitting L" w -i- 7? into linear operators 

£f : //o^(n)-//"^(n), w^-af Aw-i-fof -Vw-i-cf 

/4(n)-//-^(n), w^-a'{^w+^■vw-^e1w, 

with continuous 

A cm x C(n,R'*) X cm, a (a", c"), 

- , - = _ (5) 

A^Cm X C(0,IR'*) X C{0), a id'^,b'^,cf) 

such that c" and c? are non-negative and for some y e R and all ae A, 

r^iicriiL-+ii^riiL-. (6) 

Also find for each / a non-negative df which approximates d": df » d". These consistency conditions » Ef + /? and 
d" ~ d" are made precise as follows: 

Assumption 1. For all sequences of nodes iy^)i^\^, where in general ( = Hi) dependsoni: 
Urn supl a« - (af(yf) + 5f(yf)) + - (fof + 
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Define, for u; e J/^ (H), ^ e {1, . . . , JV = dim V;?}, 



(Ef w)/ := (yf ) (V w, V<^f > + <fof • V w + cf w, (^f >, (7a) 
(If w)^ := (yf ) <V w, V<^f > + {bf-Vw + dfw, >, (7b) 
(Q)^:=(df,<^f). (7c) 

On thie restriction to Vj- we identify tire E^w and I" w withi thieir matrix representations with respect to the nodal basis 
{<p^i}f Similarly the nodal evaluation operator corresponds then to the identity matrix Id. 

Definition 2. An operator F :V ^ is said to satisfy the Local Monotonicity Property (LMP) property if for all v £ Vi such 
that V has a non-positive minimum at the internal nodeyf, £ e {l,...,N}, we have iPv){ < 0. The operator F satisfies the 
weak Discrete Maximum Principle ( wDMP) provided that: 

if[Pw)^>:Oforall£e{\,...,N}, then minw >mm{mmw,0]. (8) 



More explicit alternative formulations of the wDMP are discussed, for example, in [5] and [6]. Note that also Id and 

satisfy this LMP property. It is clear that if P satisfies the LMP and v e Vi has a negative minimum at the internal node 
then ((F + e Id) f)^ < for all e > 0. This implies for all £ > that f + e Id satisfies the wDMP 

Assumption 2. Assume for each ae A that E" restricted to Vi has non-positive off-diagonal entries. Let hi be small enough 
so that all /z, E" - Id are monotone, i.e. so that all entries of all hi E" - Id are non-positive. Assume that for each ae A that 
I" satisfies the LMP property. 

Obtain the numerical solution Vi (T, •) e V/ by interpolation of vj- Then Vi {s^, •) e V9 at time s^ is defined, inductively, by 

-di vi{s\, •) + sup(Ef i;,(sf •) + 1? Viis^i, •) - Cf ) = 0. (9) 

a 

If all I? are then (9) is an explicit scheme, otherwise implicit. Notice that the monotonicity assumption on hi E?^ - Id is a 
time step restriction if E? has positive diagonal entries. 



Ill Well-posedness of the discrete HJB equations 

Let a^^'^^w) be a control a which maximises 

sup(E>(.f-i,.) + l>(5f,.)-q)^. (10) 
Let l^"' and E^"' be the matrices whose /th row is equal to that of 

I . ' and E. , 

respectively. Also let the /th entry of C^"' be 

Thus, informally speaking, the E.^/'" , l^"' and d^/'" are gained by 'reshuffling' the rows of the E?, I? and C", respectively. 
Notice that the maximising control in (10) may be non-unique. 

Where no ambiguity can arise we simply write and C^" without explicitly referring to k. We will make use of the 

partial ordering of IR": 

forx, yelR"; x>yif andonlyif x^- >y^, V/e{l,...,n}. 

For a collection {x"}^^^ c B?", we define the operator sup^^^ componentwise: (sup^^^x")^ = sup^^^x". The following 
lemma shows that in the linear case the wDMP turns, for functions which vanish on the boundary, into an M-matrix 
property. 
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Lemma 1. The matrices hiE^"' -\A are monotone. The matrices of hi\^"' + \A restricted to Vf are invertible diagonally 
dominant M-matrices for all w £ C([0, T] x Q). The operators andhjlf"' + Id satisfy the IMP and wDMP, respectively. 

Proof. Monotonicity of hiE'f'"' - Id is a straightforward consequence of the non-positivity of the entries of ^j E" - Id for 
all a £ j4. The LMP property of I" for the node yj only imposes a condition on the ^-th row of the matriK of I". Hence it 
is easily checked that the if"' and the hi\f"' + \d, which are composed row-wise from the I" and hi\" + Id, satisfy the LMP 
and wDMP respectively when all I? satisfy the LMP property. 

The LMP property also implies that the matrix representations of the I? restricted to V9 are weakly diagonally dominant 
for all a e A. This is because taking v = - 'Lf^i yields 

N 

^ -on..- zon/y- 
using the fact that v attains a non-positive minimum at each internal node. For j ^( the hat function attains a non- 
positive minimum at y^, giving (If ^ 0. This yields 

Therefore hi\'^''" + Id restricted to V9 is strictly diagonally dominant and thus invertible. Furthermore, since [hi\^'"' + 
Id) H- eld is similarly invertible for all e > and all off-diagonal entries are non-positive, [18, p. 114] shows that hil'^'"' + Id 
restricted to V9 is represented by an invertible M-matrix. □ 

Corollary 1. The non-linear operators if — if"" if and w — 1*^'"' + Id) w satisfy the LMP and wDMP, respectively. More- 
over, w-^ -{hiE'^'"' -\d)w is positive: ifw>0 f^en-(/z,E^"'-ld)u'> 0. 

We record a constructive proof of existence of a solution Vi e S; x V° to (9) for all fc e {0, 1,2, . . . , T/hi - 1] with the below 
Algorithm 1. This algorithm, which can be traced back to [19], is found in the continuous setting in [24] which provides 
the proof of convergence and existence of solutions. In [4] it is shown that in the discrete setting it is a semi-smooth 
Newton method that converges superlinearly. 

The algorithm to solve the non-linear problem (9) at a given time level is the following. 

Algorithm 1. Given keM and Vif^sf'*'^,-) £ V9, choose an arbitrary a £ A and find wq £ V9 such that 

[hi[f + Id) wo = hiq - [hiE'l - Id) Vi isf^\ •). 

Forme {0, 1,2,...}, inductively find Wm+i £ Vf such that 

[hill''" + Id) Wm+l = hiCf"' - [hiEf"' - Id) Vi(sf^\-). (11) 

Theorem 1. The numerical solution vi exists, is unique, solves the linear systems 

{hi\Y' + \d)Vi{s1,-)^-{hiE\'"^-\d)Vi[s^^\-) + hiQ'i' Vfc£{0,l,2,...,;i:-1}; (12) 

and is non-negative. Given k e {0,1,2,...} and 1^,(5^^"^^, •) £ V9 , the iterates of Algorithm 1 converge superlinearly to the 
unique solution y,- (s^ •) of (9) ; lim^ Wm = y,- {s^, •) . Any numerical solution vf of the linear evolution problem associated 
to a fixed a with homogeneous Dirichlet conditions, that is 



(/z,-lf -i-ld)i;?(sf,-) = -(^iE«-ld)i^f(sf+i,-) + ^,C«,fce{l,2,3,...}, with vf(^T,-) = ViiT,-), 
is an upper bound: i/,- < vf on S, x Q 



(13) 

Vi {sf,-)€V° /or a// fc £ {0, 1, . . . , I: - 1} ; 
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(a) Consistent mesh (b) Inconsistent mesh 



Figure 1: (a) illustrates a mesh that leads to a FEM discretisation of the Laplacian that is pointwise consistent with respect 
to the interpolant. This is no longer the case for the mesh depicted by (b). In (b), Ti denotes the upper-right element. 

Proof. [4, Theorem 2.1] shows existence and uniqueness of a solution Vjis^,-) given k and Vi^s'^'^^,-) and superlinear 
convergence of the algorithm: their Assumption (HI) is ensured by Lemma 1 and their Assumption (H2) is guaranteed by 
equation (5). Existence of a solution i/,- is then obtained by induction over k. 

Also Vi > on S, X n follows from induction over k. By asstmiption Vi(.T,-) > on Q. Since all entries of hiE^' - Id are 
non-positive, all entries of C^' are non-negative, and (s^"'"\ ■) > 0, (12) shows 

(hi\';' + \d)Viisf,-) = -ihiEl' - ld)i/,(5^\-) + hiCl' > 0. 

Hence by inverse positivityof hj]".' +ld, we deduce that Vi[sj,-) >Oonn. 

Finally, we prove that Vi < for all a £ A. Fix a e A. Firstly, ViiT,-] - v"{T,-]. For given A; e N assume that fK-s^^"""^,-) < 
i;f (sf+\-). Then (9) implies 

[hi\f + Id) vdsfr) < - [hiE^ - Id) vdsf^'r). 
Subtracting (13) from the above inequality and using monotonicity of hi E" - Id yields 

{hi\1 + \d)[vi(sl-)- Vf (Sf , •)) < E« - Id) ( Vf (Sf •) - Vi (5f +1, •)) < 0. 

Thus by inverse positivity of hi\f + Id we conclude that Vi (sf, •) < vf is'^, ■) on n, which completes the induction. □ 

IV Consistency properties of elliptic projections 

The Barles-Souganidis argument requires the existence of a projection operator onto the discrete function space that 
satisfies two properties. First, the projections of a smooth function must be convergent in a sufficiently strong sense, 
for example in W^'°°. Second, the discretisations of the partial differential operators must be pointwise consistent when 
applied to the projections of a smooth function, i.e. the values of the operators applied to the projections converge to the 
values of the continuous operator applied to the smooth function. In the context of classical finite difference methods, 
the interpolant to the grid satisfies these properties trivially because the operators are designed to be consistent with 
respect to interpolation. However, in the case of FEM, the nodal interpolant may fail to satisfy the consistency condition, 
even for reasonable meshes. We illustrate this behaviour in Example 1. 

Example 1. For a fixed point x in a domain, consider two sequences of meshes, such that the elements neighbouring x areas 
depicted in Figure 1. Denote (pi and (pi the -normalised hat functions associated with the node x for the meshes depicted 
respectively by (a) and (b). Let w beasmooth function; lefiaW and3bW be the nodal interpolants ofw respectively on the 
two meshes. We show that the mesh type of (a) leads to a FEM discretisation of the Laplacian that is strongly consistent with 
respect to interpolation, whereas the mesh type of(b) does not. 
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For the mesh of Figure 1 (a), it is well known that the FEM discretisation of the Laplacian coincides with a finite difference 
discretisation and that 

(V3aw, V(Pi} = — ^ {4w(x) - w(y^) - w(yf ) - w{yp - w{y^)) = -Aw{x) + 0((Ax)f ). 

(l\X} • 

For the mesh of Figure 1 (b), we sketch the calculation: first we have 

II II 2 p ^ 3 „ 1 iw{y^)-w{x)\ 

ll'^'MI 3^ h' f'^^i 2(Ax]^\-l ' {Ax]i\w(y^)-w(x) ' 



thus 



3 

VOhW- V<pi dx = (2ii/(x) - w(y?) - w{y^)) . 



Doing a similar calculation for the other elements shows that 

{VObW,S/<pi) = — ^ {4w[x) - M/(y?) - u/(yf) - u/(yf) - u/(yf)) = -|am;W + 0((Ax)?). 

^(l\X) j ^ 

We overcome this difficulty by using a different projection operator in the Barles-Souganidis argument. Given w e 
C([0, T],H^ [Q.)], denote by P i w a linear mapping into [0, T] x V,- which satisfies for all (^f e V9 

(VP,ii/(f,-),V(^f) = <Vii/(f,-),V(^f> Vfe[0,r]. (14) 

Notice that P,- coincides with the classical elliptic projection of the Laplacian if P,- w is chosen to interpolate w on the 
boundary. 

Assumption 3. There are mappings Pi satisfying (14) and there is a constant C > such that for every w e C°°(IR^) and 
i€N, 

||PiW||^yi,oo(Q) < C|| w||pvi,oo(Q) and lim IIP; w||^yi,oo(n) = 0. (15) 

i— »oo 

The conditions under which the above assumption holds for the elliptic projection typically include a condition on the 
mesh grading and on the domain. In [12], it is shown that (3) holds when O is a bounded convex polyhedral domain 
in K^, d e {2,3}, when the mesh satisfies a local quasi-uniformity condition and when the test functions vanish on the 
boundary. To apply the result for non-convex domains O and general w e C°°(IR x W^), consider for example a convex 
polyhedral domain B containing O and assume there is a locally quasi-uniform mesh on B which coincides with the 
original mesh on Q. Let be a smooth cut-off function with compact support in B such that tj s 1 on 11. Then the 
classical elliptic projection on B, acting onrjw : B U., has the required properties. Given this construction for P,, it is 
natural to refer to it as an elliptic projection. 

Lemma 2. Let w e C°°(IR x 0?'*) and let {s'^^'^°Z^ tend to t e [0, T). Then 

lim diPiW{s^^'\-)=dtw{t,-) inW^-°°{Q). (16) 

Proof. By linearity of P,- and (15) , the result follows in the limit i — oo from 

WdiPi w{sf'\-) - dtw(t, ■) \\w-^.^[a.) 
^\\diPiw(sf'\-) - diPiw{t, ■)llwi,<x,(n) + M;PiW(f,-)-Pi<3fW(f, ■)llwi,<»(n) + ll-Pi<5t ■)-dtw[t,-)\\w^'<^{n) 
<C||d;if(sf''',-)-rfiif(r,-)llivi,co(n) + C||d,-if(f, ■)-dtw(t,-)\\w^.^{Q.) + \\Pidtw[t,-)-dtw(t,-)\\w^''^(n)> 

where dfPiWit,-) = [PiW{.t+hi,-)-PiW{t,-))lhi, assuming i is sufficiently large to ensure t+hi<T. □ 

Lemmas. Lef m; e C°°(IR x IR<^) andZeflsf'''}^^ fend to f £ [0, T], {yf'''}~j tendtoxeD.. Then 

lim fEfPiM;(sf"'^"'\-) + lfP/w(5f"\-)-Cf] = L" w(t,x) - d"[x) uniformly over all a £ A. (17) 
oo V ' ' II I ;^(,-) 
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Proof. For ease of notation, the dependence of k and ( on / is made implicit. From the definition of Pi and integration by 
parts, 

I ^rCyf ) <^Pi wis^i' V0f ) + flf (yf )(VPi ii/(sf +\ •), V0f > - a" (yf ) {S!w{t, •), V(^f ) | 
= I (yf ) <V w(sf , ■), V0f > + flf (yf ) <V w(sf •), Vc^f ) - a" (yf ) (V w(r, ■), V(^f > | 

<|(a"(yf)-^f(yf)-af(yf)) <-Au;(r,-),0->| + (yf)<AM;(r,-) - Au;(sf,-),<^f>| + |af (yf)<Au;(r,-) - Aw(5f+\-),0f>| . 

Using Assumption 1 and the continuity of w: [0, T\ >— C^(Q) together with uniform boundedness of {|«"(yf)|}aeA and 
{| a?(yf) |}^^^, we conclude that 

lim sup I (yf ) (VP; u;(sf , •), V(^f > + a? (yf ) (VP,- u;(sf +S •), V(^f > - a«(yf ) (V m;(?, V(^f > | = 0. 

Owing to the Heine-Cantor theorem for all £ > 0, there is a5 > such that \lS.w{t,x) - S.w{t,y) \ < £ if |x-y| < 5. Since, for 
i sufficiently large, the support of (^f is contained in the ball B (x, 5) and since || c^f || = 1 as well as c^f > 0, we find 

|Aw(r,jc) - (Aii;(«-,-),<^f>| <£. 
As is an equi-continuous family of functions, we conclude that 

lim sup a"(yf)(AM/(f,-),<^f>-fl"(x)Ait/(f,x) =0; 

thus showing that 

lim sup I ^"(yf) (VP/ H;(sf,-),V<^f) + a?(yf)(VPi w(5f"'\-), V(^f) - (-a"(x)Au;(r,x))| = 0. (18) 

Using Assumption 3 and regularity of w, wesee thatP,M'(sf^,-) andP,- w(sf"^^,0 converge to w{t,-] in W^'°°(Q). It can then 
be shown by analogous estimates and by using the equi-continuity of {c"},^^^ and that 

lim sup I {hf ■ VPi u/(sf , •), (^f ) + {bf ■ VPi ii/(sf +\ O.f^f > - b"{x) • Vit/(f, x) I = 0, (19a) 

lim sup I (dfPi w{sf, ■),^\) + {cfPt M;(sf+\ O.f^f > - c"(x) w{t,x)\ = 0, (19b) 

limsup|(d?,<^f>-d"(x)Uo. (19c) 

Combining equations (18) and (19) yields (17). □ 



V Sub- and supersolution 

Set 

v*{t,x)- sup limsup t;/(Spyf), v*(t,x)- inf liminfy, (x^yf) 

where the limit superior and limit inferior are taken over all sequences of nodes in [0, T] x n which converge to (t,x) e 
[0, r] X Q. By construction, v* is upper and i"* lower semi-continuous. With the use of elliptic projection operators key 
steps of the convergence proof in [3] , which is stated there in a suitable form for finite difference methods, are transferred 
to finite element schemes, which do not satisfy the consistency condition in [3] . 

Theorem 2. The function v* is a viscosity subsolution of (4) and is a viscosity supersolution of (4) . 

Proof. Step 1 (v* is a subsolution). To show that v* is a viscosity subsolution, suppose that w e C°°(IR x IR'*) is a test 
function such that v* -w has a strict local maximum at ( s, y) e (0, T) x Q, with v* is,y) = w[s,y]. Consider a neighbourhood 
B:={it, x]£iO, T]xD, : \t-s\ + \x-y\<S} with S>0 such that 

v*is,y)-wis,y) > v* {t,x) - w(t,x) V(f,x) eB\(s,y). 
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Choose i sufficiently large for B to contain nodes. Let isf.yf) denote the position where Vi[s^,yf) - Pi w{s^,yf) attains 
the maximum among all nodes £ B. Let us pass to a subsequence {(5f(j).yfy))}yeM of yf^l/eN which 

{i';(^fy).yf(y))}yeM Converges to the limit superior of {vi(.sf,y^^)}i^^. By compactness of B, there is a subsequence of 
{(•sf(y).yf(j))}yeM Converging to a point (5,y) e B. Then ^/"^(^fyj-yfyj) wis,y) from (15) and by continuity of w. As the 
(Sj^, yj ] are maximisers, one has 

v*(.s,y) - w[s,y) = limsup i//(sfyj,yfyj) -P,-M;(sfyj,yfyj) = v*(s,y) - w(s,yy, 
i— »oo 

hence is,y) = (5,y) since (s,y) is a strict maximiser of v* - w on B. Thus there is a subsequence of maximising nodes 
converging to is, y) to which we now pass without change of notation: (sf, yf ) — is, y) . It follows that 

Vi isf, yf ) - Pi wis^, yf)^v* is, y) - wis, y) = 0. (20) 

Moreover, because of (sf.yf) (s,y), the neighbours of the (s'^,yp eventually also belong to B: For / sufficiently large we 
have (s^, y^) eiJif?c£{fc, A;+l} and yj e supp (^f ; in which case 

Viis'^,yf]-PiWis'^,yf) < Viis'^,yf)-Piwisf,yf) o Piwis'^,yf)+Hi > J//(s^,yf), 

with Hi = ViiSf,yp-PiwiSf,yp. Notice that m — as / — oo because of (20). 

Since the matrices E?^ have non-zero off diagonal entries (E?)^^ only if y^ e supp(^f , we have for all a e A 

(ihiEf - Id) [PiWis''^\-)+i,i])^ < (ihiEf - \d)Viisf^\-))^. 
By the LMP property and linearity of I", since P,- wisf, •) + m - Vi is^, •) has a non-positive minimum at yf , 

[ihi If + Id) \Pi wis'l, •) + Mi] )^ ^ ((^/ 1" + Id) Vi (sf , ■)]^. 
From the definition of the scheme, 

= -rf/i;,(sf,yf) -HsupfEf i;,(sf+i,-) + l"i^K5f,-) -Cf) . (21a) 

> -di [Pi wis\, yf ) + ju,) -I- sup (e? [Pi u/(sf +1, ■) + Hi] + \f [Pt w^Csf , ■) -i- m/) - Cf )^ (21b) 

= -diPi u/(5f ,yf ) + sup [(e«P,- w{s\^\-] + l^Pi wisf, ■) - C«] + ^i; (cf + 5f )1 (21c) 

> -diPi w{s\, yf ) + sup (e?P/ w{s\^\ ■) + I'fPi wisf, •) - C?) - y 1^,- 1 . (21d) 



Since 



sup (Ef P; wisf+h, •) + \fPi w(sf, •) - ) - sup (L" M/(s, y) - ci«(y)) 
< sup I (e^'P; u;(sf •) + \^Pi wisf, •) - ) - (I" wis, y) - d^(y)) I , 
Lemmas 2 and 3 show that after taking the limit / — co in inequality (2 Id) and recalling that jU, — 0, we obtain 



> -dt wis, y) -I- sup [L" wis, y) - d" (y)) . (22) 



Therefore v* is a viscosity subsolution. 



Step 2 (v* is a super solution). Arguments similar to those above show that i/* is a viscosity supersolution, where the 

principal changes are that one considers w e C°°(IR x U'^) such that Vt, - w has a strict local minimum at some (s,y) e 
(0, T) X n with V* is, y) = wis, y) . Using analogous notation, inequality (21d) corresponds to 

< -diPi wisf, yf ) + sup (e?P/ wisf^\ •) + l?P/ wisf, ■)-Cf] +r\ni\, 

i.e. there is a slight asymmetry in the argument due to the last sign in (2 Id). Nevertheless it is then deduced that 

< -5f m;(s, y) -I- sup (i" wis, y) - (y)) . 

Thus I'* is a viscosity supersolution. □ 
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VI Uniform convergence 



We now turn to the initial and boundary conditions. Together with the sub- and supersolution property we appeal to a 
comparison principle to obtain uniform convergence of the numerical solutions. 



For each aeA, define 



j;"'*(r,jc)= sup limsupffCsf.yf); 



where the vf are as in (13) and the limit superior is taken over all sequences of nodes which converge to (f, x) e [0, T] x Q. 
Assumption 4. Suppose that for each (f, x) e [0, T] x dn 

inf i;"'*(?,x)=0. (23) 

Before further considerations, let us motivate Assumption 4 with a simple example. As a side remark, this example also 
illustrates how in some settings Kushner-Dupuis finite difference schemes, as described in [23, 17], maybe interpreted as 
finite element methods in the framework of this paper. 

Example 2. Consider the backward time-dependent equation in one spatial dimension 

-Vt + \Vx\ = l on(0,l)x(-l,l), (24) 
with boundary conditions v = on [0,1] x {-l, 1} u {1} x [-1, l]. Equation (24) may be re-written inHJBformas 

-vt+ sup {aux-V)-Q. 

ae{-l,l} 

The viscosity solution is v = min(l -t,l-\x\). We choose a uniform mesh with element size2(Ax)i and we use a fully explicit 
discretisation, where monotonicity will be achieved by using the method of artificial diffusion, as described in [5]. Thus we 
have 

[^"w)^ = e{dxW,dx(p\) + a{dxW,(p\), 

where e is the artificial diffusion parameter to be chosen to obtain a monotone scheme. Calculating the entries shows that 
the E? are of the form 



f_2 i_ 

(Ax)? (Ax)? 



1 2 I_ 

■(Ax)? (Ax)? (Ax)? 



1 2 
(Ax)? (Ax)? ) 



1 



2(Ax)j 



1 



1 



2(Ax); 



2(Ax)i 



2(Ax)i 



For monotonicity we require that all off-diagonal terms of the E? be non-positive, i.e. we require e > (Ax)/ 12. For example 
the special choice e = (Ax) , /2 yields 



(Ax); 



EH 



1 



1 



(Ax)i (Ax)j 



1 1 
■(Ax)i (Ax),-^ 



L_ 

(Ax)i (Ax)i 



1 



1 



(Ax)i (Ax)( 



(Ax)i / 



This is equivalent to discretising the spatial part of- Vt+Vx with backward finite differences and discretising the spatial part 
of-V(- Vx with forward finite differences, as can be done in applying a Kushner-Dupuis scheme. It can then be deduced, 
whilst using appropriate time steps, that v] approximates the solution of 



-vt+vx^lon (0, 1) X (-1, 1), j; = on (0, T) x {-1} u {1} x (-1, 1); 
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while V. approximates the solution of 

-Vt-Vx = lon (0, 1) X (-1, 1), v^O on (0, T) x {1} u {1} x (-1, 1). 
Consequently, Assumption 4 is enforced by v^'* on [0, 1] x {-1} and by v~^'* on [0, 1] x {l}. 

Recall from Theorem 1 that 

0<Vi< v" for all a; e A, 

and note that by construction < i/* < i/*. Assumption 4 thus implies that i'*|[o,T]xan = i^*l[o,T]xan = 0. Observe that 
because (23) holds in particular for all (f,x) e {T} x 50, Assumption 4 implicidy enforces that the initial condition vt 
vanishes on dD. as the interpolate vj at the final time. 

Lemma 4. The sub and super-solutions v* and v* satisfy 

v*{T, •) = v*(T,-) = Vt on n. (25) 



Proof. Flk £ > and choose a v^e C°°(IR'') such that vj -2e > > vj - 3e. Owing to Assumption 3 there is « e N such 
that IIP; Vj - Vj\\L''°{n) ^ £ and HJ/fr - frlko^cn) s e for all / > n. Hence, for / > n, 

Vi (r, •) = J/ Vt > Pi v'j>Vt- 4£. (26) 

Recalling (14) and as v^ e C°°(IR'*), it is clear that there exists K = K(^e) > which bounds 

m + \-)PiV^-q),\^\-[c^iyi) + a^iy^)){Av^,<P^) + a^ 

for all J e N, ^ £ {l,...,iV} and a £ A Define Wi - PiVj-K{T- t). To show that viis^,-) ^ u)i(s^,-) assume viis^^^,-) ^ 
w^•(sf"^^ noting (26) for 5^+^ = T. Fix an / and ( and let a = ivi) as for (10). From 

-diWi^s\,/i) + [£^Wi{s\^\-) + \'lwi{s\^\-))l^-K+[[E'f + \'^)PiV^^^ 

< [Cf), 'L^' -divdslyf) + (Ef vdsr\.) + If vdsf'\.)), 

we may deduce that 

((/I,-!,!''- + Id) - > ((/j^E;'- - Id) - > a 

Note that Vjis'^,-) e vanishes on 5Q and ^/(sj^, ■) < on 5Q. Thus Lemma 1 and (8) imply I'K-sf. ^ i4'/(sf, ■) on Q. 
Because i<r is independent of / and Pi Vj Vj as ; oo, we have for any sequence (5f,yf)-(r,x),x£a 

liminf I// (sf,yf) >liminf u;,(spyf) > i'r(Jc) -4e. 

So i/*(r,-) > i/r -4e. Since £ was arbitrary, i/*(r,-) > i/j. The argument for showing that v* < Vt is analogous with 
Wi = PiVj + K{T-t) and i/r + 2£ < ff. < fr + 3£. To conclude, Vt < v* {T, ■) <v* {!,■)< Vt, which proves (25). □ 

The proof of Lemma 4 is related to the arguments in [17, p. 335]. In the next assumption we draw upon one of the 
building blocks of the theory of viscosity solutions, namely the extension of classical comparison principles to spaces of 
semi-continuous functions, cf. [10, Sec. 5] and [17, p. 219]. 

Assumption 5. LetD be a lower semi-continuous supersolution withv\[o,T]>cdn = andviT,-) - vt. Similarly, letv bean 
upper semi-continuous subsolution with v\io,T]xa = andv{T,-) = vt- Thenv<v. 



Let t = dsf + (,l-d)s^'^^ £ [s^^,s^"^^] lie between two time steps, d e [0,1]. Then we interpret i/((f,-) as the linear interpolant 
between Viisf,-) and vtis^^^,-): 

Vi(.t,-)=dVi(.s^,-) + a-'d)Vi{sf*\-). (27) 
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Theorems. Onehasv* = v* - v, where vis the unique viscosity solution of equation [A] withv{T,-) - vj andv\io^T]xdn = 0- 
Furthermore 

lim II Vi - J^||L«'((o,r)xn) = 0. (28) 



Proof. The previous assumption implies that v* > v* thus v* = v* = v. Select for each / e N a point (t,-,x/) e [0, T] x n 
such that 

\\vi - i'||L°°((o,r)xn) = - v\{ti,Xi). 

Such (ti,Xi) exist as Vi - f is a continuous function on a compact domain. Let x,- belong to (the closure of) the element T 
of the finite element mesh and t e[s^, s^"^^]; then viiti, x/) is a weighted average of the values of Vi at the corners of the 
slab [Sp s]^"*"^] X T. Thus there is a corner (s^, yf") of the slab such that 

\\Vi - i/||L°o((o,r)xn) < jfK^f.yf) - v(ti,Xi)\. 



If (28) was wrong we could select a subsequence and an £ > such that 

k e 



liminf I i^,-o-)(sL,yfy)) - l/(?/y),X,-y))| > £. 



7— >oo 

By possibly passing to a further subsequence we may assume that {(r, (J), x,y))}; converges to an (r,x) e [0, T] xQ. However, 
this contradicts 

v{t,x) = i',(r,x) < liminf i',(;)(sL),yfy)) < lim sup i'/(j)(sfy),yf( -J ) < v*{t,x) = i^(r,x). 
Thus (28) holds. □ 



VII Gradient convergence 

For shorthand, let W - 14^^'°° ((0, T) x Q). It is convenient to introduce the discrete spaces 

Wi := [vECao, T],vP) : v\^^k^^k*i]^n is affine in time}, 
which means that functions in Wi have between two time-steps the form of (27). Observe that Wi ciy for all j e N. 

Fix an arbitrary aeA.lt is convenient to view E? and I" as bilinear forms on (Q) x V/. Functions m e V/ have the nodal 
representation 

u(y]^J^u(yi)<pi{y). 



To test with functions other than (pj we introduce the following bilinear form as a partially discrete pivot: for we. (Q) 
and ueVi 

«E^ w, M)) := X "^yf ) («r (yf ) <^ ^'Pi^ + {bf-Vw + cfw,<p^i}). 
e 

We use corresponding interpretation for u/, u)) and also 

«M/, M» = «ldM/, M» = 2^ u/(yf) u(yf)||(/)f llinn) and «C?, u)) = uiyf) {d%<l)\) = (df, u). 

Assume that for the chosen a: 

lTlhi)-\ 

I "^^iHlo.Tim im ~ [ii'^i - Id) iy(sf +\ ■) + (/i,- If + Id) w{sf, ■), M^(5f , •)))) + i « t^(r, •), w(r, •)» + II iv{T, ■) ii^ij^j 

[Tlhi)-l 

= L k•(<Efu/(s^^o+lfM/(s^o,M/(s^o))+i<(M/(5f+^o-u;(5f,o,Msf+^o-u/(s^^ (29) 

+ |«u/(0,-),it/(0,-))) + ||M;(r,-)ll^i,^) 
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for all M/ e Wi with w>0 and i e N, where (*) is a simple reformulation in terms of a telescope sum. 
Due to the definition of the numerical method and the non-negativity of the Vi, 



k=0 
iTlhi)-l 



< ^ {{hiq, Vi(.sf. •))) + i {{vdZ ■), vdT, ■)}} + II vdT, ■) ll^ij^j 



fc=0 

< rildf llinn, ||i;Hli-ao,rixn) + lk,-(r,-)ll^i(n,. 

Thus, with the L°° control established in the previous section, it is apparent that the Vi are bounded in L^([0, T],H^ (O)) 
provided that vi(.T,-) = Ji^r are bounded in H^{D); this condition holds if v{T,-) e iy^'°°(n). The first convergence re- 
sult for the gradient is therefore that, owing to the Banach-Alaoglu theorem, vi v weakly in L^{[0, T],H^{Q]), using 
L°°((0, T) X D) convergence to pass from L^([0, T],H^{P)) weak convergence of subsequences to L^([0, r],//^(n)) weak 
convergence of the whole sequence. 

The question arises under which circumstances the convergence in the gradient is also strong. We demonstrate this 
under the below Assumption 6. We note that supposing (29) points towards uniform ellipticity of L". Let Aq be the level 
set { ( f , x) e (0, T) X Q : ( x) = 0} . For a smooth v the boundary of Aq is always ad -I dimensional set if is a regular 
value. 

Assumption 6. The value function v belongs to the space W = W^'°°{{0, T) x Q) and the d -dimensional Lebesgue measure 
of the boundary of Aq vanishes: volidAo) =0. The coefficients df^ andd" belong to W^'°°{n.) and [29) is satisfied. 

Let us suppose momentarily that there are approximations Q, VEWfto v such that Q,- v < i/, for all / e N and 

lim \\v-Qiv\\L2^o,T],HHn))^^' 

and 

(Tlhi)-l 

lim X {{ihiE'}-\d)QiV{sf^\-] + {hi\'^ + \d)QMsl-]Avi-Qiv)isfr)}-0. (30) 
We will construct such Q, v below. With = i;,- (5^, ■) - Q,- v{s^, ■), 

[Tlhi)-l 
k=0 

(.Tlhi)-1 {Tlhi)-\ 

= E i[hiE'}-\d)vi{sf^\-) + ihi\'f + \d)vi{s1,-U''}- L {{ihiE'^-\d)QiV{sf^\-] + ihi\'f + \d)QiV{sf,-U'} 

k=0 k=0 

W 2: ihiq,^'}- L {{{htE^-\d)QiVis^^\-) + {hi\f + \d)QiV{sl-)A (31) 

k=0 k=0 

using in {*) the ntmierical scheme, f = and that, due to the assumptions on the Q,, the sign of Vi - Qi v is known. 
Since 

{Tlhi)-1 i.Tlhi)-l 

ihiq,e))^\\df\h2^n) L hi{}\vds\r)-v{s\MLHn) + \\v{s''tr)-QiV{sl-)\\L2^a)) 

k=0 k=0 

^WdfhHn) i\\Vi-v\\L2mT)xa) + \\v-Qii^\\LH{o,T)xn))' 

the first term in (31) vanishes as / — 00. The second term vanishes due to (30). Hence \ vi - i^li2([o,ri,Hi(n)) — as i — 00. 

Theorem 4. If there is anaeA such that Assumption 6 holds, then the numerical solutions converge to the exact solution 
strongly in ( [0, T] , (Q)) . 
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Proof. It remains to show that suitable Q; can be constructed, given Assumption 6. Denoting the nodal interpolant on 
[0, T] X n by J, we define 

Qi: W-^Wi, w-^3imsy.{w-\\v-Vi\\L-onfi,T)^n),0]. (32) 

Observe that the max operator in (32) switches between the first and second argument in the vicinity of (5Ao for ; suffi- 
ciently large. Furthermore, Q/ v £Wi satisfies homogeneous boundary conditions and QiV < y,- and, by the mean value 
theorem, 

IIQjf llwi.oo((o,r)xn) - ll*'llwi''»((o,jnxn)- 
Note also that for all nodes and time levels 

< {Vi - QiV) (sf ,yf) = min{(i;/ - v) (sf ,yf) + II vi - I'lli^ao.nxn) - MSi.yf)} < 2 1| vi - y|lL<x.((o,T)xn) • 

Consider the set T; of points which is not 'affected by the cut-off below 0' in (32) in the sense that 

r, := Ut,x) e (0, T) X n : infQ/i;(?,x) > or {t,x) e Aq}. 

The set r^. contains the points which are at least one element's length away from the boundary of \ 5Ao: 
r;. := {(r, X) e r,- : Us, y) e (0, T) x n : || (r, x) - [s, y) || < sup hj + iAx)j] e r,- \ dAo}. 

Notice that r,- andr^. are hierarchical families. Since ||i;- Vi\\Loo[^[Q^T)xQ.) and hi + (Ax), — as / — oo it follows that 

Ur; = ((o,r)xn)\aAo. 

Crucially, {dtQjV)\Y'. = (.diOjV)\Y'. and {VQjV)\p. = (VJyi/)lr'. for; > i. 

For every e > there are i, j e N such that vol(0 \ Fp < and II Qjt f - II //i (r' ) < e for all fc > / . Therefore 

< II Qfci^ - I'll HHT'i) + \/vol(n \ T'.) II v\\ 

vvi.°°((o,r)xn) ^ + lli'llwi'«'((o,r)xn))> 

giving strong convergence in //^((O, T) x Q), meaning convergence in the spatial gradient and the time derivative. The 
terms connected to the time derivative in (30) vanish in the limit as 

iTlhi)-l [Tlhi)-\ 

L iQiV{s\^\-)-QiV{s\,-),^^)) = i: hii{dtQiV)\ ,^,,,ye))<\\dtv\\i2(^^,T).a) U''\\LHif>,T).ay 

Recall that 

<<lfQ/l'(5^0,<*^>> = Z(l'^-Q^l')(sf,yf)(^^'(yf)<VQ,■l;(sf,0,V(/.f) + <Sf•VQ,•y(5f,0 + a« Q,i'(sf,-),0f». 
e 

The lower-order terms vanish due to the uniform convergence of i/,- - Q, f to and the bound 

sup II • VQ,- i/(sf , ■) + ef Qi i;(5f , •) llioo(n) < oo. 
i 

We note for the second-order term that 

^(i/,-Q;i;)(sf,yf)af(yf)(VQ,i/(5f,0,V(^f)= <VQ;i;(sf,0,VJ,(^r(i',-Q/i'))(5f,-)>, 

e 

so that in (30) the implicit part of the second-order term becomes 

X; hi{VQiv(sf,-),^Oiiaf{Vi-QiV)){sf,-)}= / {3i^QiV,3i^0iaf[vi-Qiv)))dt, (33) 



14 



where 3j maps any w : [0, T] L^(0;IR'') onto the step function with (3,- w)\^^k ^t+ij = w[s'f,-). Note that 3; VQ,y converges 
strongly in ((0, T) x Q; R'*) . At a time g [0, T) the bound 

IIVJK^fd'/ -Qif^))lli2(n;Rd) < IIVJ,-(^f j^;)||i2(n;Rd) + ll^f Qit/|lvvi,<x,(n) < U°'\\w^.oo^n) ■ (II t^ilbi(n) + lli^llwi.~(n)) 
follows from an inverse estimate and 

The convergence 

lim r (w,diV3i{af{Vi-QiV)))dt=-]im {V ■ w,3i0iianvi-QiV)))dt = O 

oo Jo oo Jo 

with test functions w in the dense subset C^((0, T) x 0;IR'') gives weak convergence of VJiid^iVi - QiV]) in L^((0, T) x 
Q;IR'*), see [28, p. 121]. Combing weak and strong convergence [29, Prop. 21.23], it is ensured that (33) converges to as 
j — CO. A similar argument guarantees that /z,- «E?Q,- v^.sf'*'^,-),^''}} vanishes in the limit. □ 

The regularity of the exact value function v is, for instance, discussed in Section IV8 and IV.9 of [17]. Another item of 
Assumption 6, namely the justification of (29), is examined in the following example: 

Example 3. a) Suppose that a" is positive and constant and, for all smooth w, 

L"w = l"w=-a"'Aw + b"'-Vw + c"w, E"w^O, 
and, to obtain semi-definiteness in the lower-order terms, c" - 1 V • fo" > 0. Then, for weWi, 

(Tlhi)-1 (Tlhi)-\ 

«"l"^li^(loriHi(n))^«" L hi{S7wisf,-),Vw(sf,-)) + \\wiT,-]\\l,^^^^ X hi{{\'}w{s1,-),w{sl-]} + \\w[T,-nl,^^y 
k=0 k=0 

b) Suppose that a" e W^'°° (O) is non-constant, positive, uniformly bounded from below and that - | (V • fo"^ + Aa") > 0, 
noting for smooth w: 

(L" w,w) = {a"Vw,Vw} + {{c"-^{V-b" + Aa"))w,w}. 
Again choosing a fully implicit scheme with L" - I", the highest order term in {{\"w, w)) is at time s'^: 

X w[sl yf ) (sf , yf ) <V wisf, •), V</)f ) = (V wisf, •), VJ,- (a" (sf , •) wisf, •))). 
e 

According to Theorem 2.1 in [11] there is a constant C = C(||a"||^v2,oo(n)) such that fori sufficiently large 

{\Iw,\I3i(a"w)) - (Vw.Va^w) < W'^wWj^Zf^nM'') ' H^K^""^) - ^""'llffHn) ^ C (Ax); II w|l^i(0)' 

using that the rj appearing in the proof in [11] is defined in terms of nodal interpolation. It then follows from Poincar&s 
inequality that there is some C such that for C {Ax) i < iinfna" thatlwl^^^^^.^ < {(\"w, w)) for w e Vf, implying (29). 

VIII Example: the method of artificial diffusion 

The purpose of this section is to provide a way of constructing the operators E" and I" in order to satisfy Assumptions 1 
and 2. This approach, called the method of artificial diffusion, is based on the fact that for stricdy acute meshes, the dis- 
crete Laplacian is monotone. Further details on the method of artificial diffusion and monotone finite element schemes 

may, for example, be found in [5], [9] and [27]. 

Let 7 1 be the mesh corresponding to the finite element space V;. Given a function / : Q ^ IR'* we denote 

l/lr:=(LI|/j|l™)'' r£T,-,/£N. 

7=1 
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If / is elementwise constant then |/| j- is simply tlie Euclidean norm of / on T. Let {IS.x)t denote the diameter of T. We 
assume that the meshes T, are strictly acute [5] in the sense that there exists d e i0,nl2) such that 

V^f •V(^||j,<-sin(^)|V</)f|r|V(/)?|r V/,/<iVVieN. (34) 

We choose a splitting of the form a" = a" + a", b" = + b", c" = c" + c" and d" = d", where all terms are in C(0), a" 
and are non-negative and all c? and c? a 
such that for all T which have as vertex: 



and 5? are non-negative and all c? and c? are non-negative and satisfy inequality (6). Choose non-negative v"'^ and v"'^ 



(I fo? I r + (Ax) r II c? || i<-(r) ) < vf sin(fi) | Vc^f | r voKD , (35a) 
(|Sf|r+(Ax)rll^?llL-{r))<^f'^sin(-9)|V<^flrvol(r). (35b) 

Choose df and df both in C(n.) such that a"(yf) > max{a?(yf),v"'^} and ^"(yf) > max{^?(yf),'9"'^}. Now suppose that 
If e V/ has a non-positive minimum at an interior node yf. By extending the arguments of [5] , we show that 

{Efw)f<0, [\fw)f<0. (36) 

We illustrate the proof of (36) for the implicit term. From the strict acuteness condition on the mesh, it can be shown that 
on the restriction to T [5, Lemma 3.1] 

Vw-V(p\ = cos(Z(Vw,V(/)f))|Vu;|r|V^flr<-sin(i»)|VM;|r|V(/)flr. 

Using a? > 0, wiy^) < and ll^f IliU^) = 1, 

{5fw,<p^i)= [ 5f(x)(u/(yf)+Vii/(x)-(x-yf))(^f(x)dx 
Jn 

< [ 5f(x)VM;(x)-(x-yf)(^f(x)dx<Xll^rilL-(r)|Vw|T(Ax)r. 
Ja T 

Consequendy, 

(.\fw)e = df[y^i){Vw,\/(p<} + (bf-Vw + dfw,<p^i) 

^ L -^f fJ'f ^ sin(^) 1^ "^1 r I V<^f I r voKD + 1 I r I V wl r + II lli-m I V w| r (Ax) r 

T 

<ZlVw|r((l^flr+(Ax)r||ef||ioo(r))-'^f^ sin(-9) |V(^f IrvoKD) < 0. 

T 

The proof of (E" w)^ < is analogous. As hat functions (/if attain a non-positive minimum at all yj where j (, all off- 
diagonal entries of E" are non-positive. Hence with a suitable time step restriction the ^;E? - Id are monotone, which 
ensures that Assumption 2 is satisfied. 

The scaling of the terms in (35) with respect to (Ax) r leads to Assumption 1. Due to shape-regularity all elements T on a 
patch are of comparable size; giving ||^f llmnj ^ Cvol(r) for all r c supp</)f with a constant C which is independent of h 
and (. Hence in (35), we see that 

^ voi(r) ^ 1 



vol(r)|V^^lr> 



(Ax)rl|(/)fllii(n) C(Ax)t 



Thus, if v"'^ and v"'^ are chosen optimally then for T c suppc^f 

v"'^ = 0(sup{|fc^'lr(Ax)r + llcf||i«=(r)(Ax)2j.}), = 0(sup{|Sf |r(Ax)r + ll^f ||L~(r)(Ax)|}). (37) 

T T 

With (37) in mind we return to the time step restriction for semi-implicit and explicit methods. The non-positivity of the 
diagonal terms of hi E" - Id expands to 

l>hi(a'f{yi){V<pi,V^l) + {b'f-v4 + c't4,^l)] = hi[^^^^ 

Therefore the time step restriction imposed by L" is hi < supj'((Ax)y/a?(yf)), yf e T, if there is a non-zero df and i is 
large. It is hi < supj-((Ax)r/|fo?(yf)lr) if all df = 0, / e N, and there are non-zero b", and is 0(1) if all df and b" vanish. 
There is no restriction if also all c? are zero. 
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